IODS course project
Anna-Liina Olkinuora
About the project
Tällä kurssilla käytetään ohjelmistoja Rstudio, Github ja Github desktop. Rstudiota käytetään tilastolliseen analyysiin ja Githubia tiedostojen säilömiseen. Github desktopilla yhdistetään tietokoneen tiedostot ja säilytyspaikka Githubissa.
https://github.com/aolkinuo/IODS-project
Regression and model validation
learning2014=read.table("data/learning2014.txt", sep = ",", header = TRUE)
str(learning2014)## 'data.frame': 166 obs. of 7 variables:
## $ gender : Factor w/ 2 levels "F","M": 1 2 1 2 2 1 2 1 2 1 ...
## $ Age : int 53 55 49 53 49 38 50 37 37 42 ...
## $ Attitude: int 37 31 25 35 37 38 35 29 38 21 ...
## $ deep : num 3.58 2.92 3.5 3.5 3.67 ...
## $ stra : num 3.38 2.75 3.62 3.12 3.62 ...
## $ surf : num 2.58 3.17 2.25 2.25 2.83 ...
## $ Points : int 25 12 24 10 22 21 21 31 24 26 ...dim(learning2014)## [1] 166 7Aineistossa on 166 havaintoyksikköä ja seitsemän muuttujaa. Muuttujat Age, Attitude, Points ovat numeerisia kokonaislukumuuttujia ja ne kertovat henkilön iän, asenteen tilastotiedettä kohtaan ja koepisteet. Ikä on laskettu syntymävuoden avulla. Muuttujat deep, stra ja surf ovat numeerisia muuttujia. Ne kertovat henkilön vastaukset syvällisiin, strategisiin ja pinnallisiin kysymyksiin. gender on luokiteltu muuttuja ja sen arvo 1 kuvaa naista ja arvo 2 miestä.
summary(learning2014)## gender Age Attitude deep stra
## F:110 Min. :17.00 Min. :14.00 Min. :1.583 Min. :1.250
## M: 56 1st Qu.:21.00 1st Qu.:26.00 1st Qu.:3.333 1st Qu.:2.625
## Median :22.00 Median :32.00 Median :3.667 Median :3.188
## Mean :25.51 Mean :31.43 Mean :3.680 Mean :3.121
## 3rd Qu.:27.00 3rd Qu.:37.00 3rd Qu.:4.083 3rd Qu.:3.625
## Max. :55.00 Max. :50.00 Max. :4.917 Max. :5.000
## surf Points
## Min. :1.583 Min. : 7.00
## 1st Qu.:2.417 1st Qu.:19.00
## Median :2.833 Median :23.00
## Mean :2.787 Mean :22.72
## 3rd Qu.:3.167 3rd Qu.:27.75
## Max. :4.333 Max. :33.00Aineiston yhteenvedosta, joka on tämän tekstin yläpuolella nähdään muuttujien Age, Attitude, deep, stra, surf ja Points minimit, maksimit, ensimmäinen ja kolmas kvantiili, mediaani ja keskiarvo. Muuttujasta gender nähdään sen kahden arvon frekvenssit.
library(ggplot2)
library(GGally)
pairs(learning2014[-1], col = learning2014$gender)
p=ggpairs(learning2014, mapping = aes(col = gender, alpha = 0.3), lower = list(combo = wrap("facethist", bins = 20)))
p
Kuvasta nähdään että muuttujat Attitude, deep ja stra ovat jakautuneet lähes normaalisti naisten ja miesten joukossa. Myös surf on naisten osalta jakautunut lähes normaalisti. Age on vinoutunut, skewed, oikealle ja myös surf on vinoutunut oikealle miesten osalta. Points-muuttujan jakauma on vinoutunut vasemmalle. Korrelaatiokertoimista nähdään että muuttujat Attitude ja Age, deep ja Age, deep ja Attitude, stra ja Attitude, Points ja Attitude, stra ja deep, Points ja deep korreloivat suhteellisen paljon naisten joukossa. Muuttujat Attitude ja Age, deep ja Age, stra ja Age, surf ja Attitude, Points ja Attitude, surf ja deep, Points ja deep korreloivat suhteellisen paljon miesten joukossa. Suuren korrelaatiokertoimen tulkinta on esimerkiksi muuttujien Attitude ja Age kohdalla se että asenteen tilastotiedettä kohtaan ja iän välillä on riippuvuutta. Jos korrelaatiokerroin olisi yksi, toisen muuttujan arvon voisi laskea tarkasti, lineaarisesti toisen arvosta.
malli1=lm(Points~Age+Attitude+deep, data=learning2014)
summary(malli1)##
## Call:
## lm(formula = Points ~ Age + Attitude + deep, data = learning2014)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.206 -3.430 0.204 3.979 10.950
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.60773 3.38966 4.605 8.32e-06 ***
## Age -0.07716 0.05322 -1.450 0.149
## Attitude 0.35941 0.05696 6.310 2.56e-09 ***
## deep -0.60275 0.75031 -0.803 0.423
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.307 on 162 degrees of freedom
## Multiple R-squared: 0.2043, Adjusted R-squared: 0.1896
## F-statistic: 13.87 on 3 and 162 DF, p-value: 4.305e-08Kun estimoidaan pienimmän neliösumman menetelmällä malli jossa ikä, asenne tilastotiedettä kohtaan ja syvälliset kysymykset selittävät koepisteitä, saadaan tulokseksi tulostus joka on tämän tekstin yläpuolella.
Ensimmäisessä sarakkeessa jossa on numeroita näkyy selittäviin muuttujiin liittyvien parametrien eli parametrivektorin arvot. Niistä nähdään että kun ikä kasvaa yhdellä, koepisteet laskevat noin -0.08 verran. Kun asenne tilastotiedettä kohtaan nousee yhdellä, koepisteet kasvavat noin 0.36 verran. Kun taas syvälliset kysymykset nousevat yhdellä, koepisteet laskevat 0.6 verran. Kun ikä, asenne ja syvälliset kysymykset ovat nolla, koepisteet saavan vakiotermin mukaisen arvon eli noin 15.61.
Tähdet kuvassa kertovat selittävien muuttujien tilastollisesta merkitsevyydestä. Tähtien taustalla on tilastollinen testi jossa testataan sitä voidaanko nollahypoteesi siitä että selittävään muuttujaan liittyvä parametri on nolla hylätä. Nollahypoteesit että vakiotermin ja asennemuuttujan parametrit ovat nolla voidaan hylätä 0.1 prosentin merkitsevyystasolla eli asennemuuttujan ja vakiotermin parametrit eli vaikutukset selitettävään koepistemuuttujaan ovat merkitseviä.Nollahypoteeseja siitä että asenne- ja syvälliset kysymykset-muuttujien parametrit ovat nolla, ei voida hylätä. Asenne- ja syvälliset kysymykset -muuttujien vaikutukset koepistemuuttujaan eivät siis ole merkitseviä. Poistetaan muuttujien vaikutusten huonon merkitsevyyden puolesta asenne- ja syvälliset #kysymykset-muuttujat mallista ja estimoidaan pienimmän neliösumman #menetelmällä malli jossa on yksi selittävä muuttuja. Se on asenne tilastotiedettä kohtaan.
malli2=lm(Points~Attitude, data=learning2014)
summary(malli2)##
## Call:
## lm(formula = Points ~ Attitude, data = learning2014)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.9763 -3.2119 0.4339 4.1534 10.6645
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.63715 1.83035 6.358 1.95e-09 ***
## Attitude 0.35255 0.05674 6.214 4.12e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.32 on 164 degrees of freedom
## Multiple R-squared: 0.1906, Adjusted R-squared: 0.1856
## F-statistic: 38.61 on 1 and 164 DF, p-value: 4.119e-09Tulostuksesta nähdään että asennemuuttujan parametrin arvo on noin 0.35. Se tarkoittaa että kun asennemuuttuja kasvaa yhdellä, koepisteet kasvavat noin 0.35 verran. Asennemuuttujaan ja vakiotermiin liittyvät parametrit ovat tilastollisesti merkitseviä 0.1 prosentin merkitsevyystasoilla.
Multiple R-squared ja Adjusted R-squared ovat mallin selitysasteita. Selitysaste kertoo sen kuinka paljon malli, eli selittävät muuttujat jotka on mallinnettu käyttämällä lineaarista regressiota, selittää selitettävän muuttujan vaihtelusta. Yläpuolella olevan mallin selitysaste kertoo siis sen että asenteen tilastotiedettä kohtaan ja vakiotermin vaihtelut selittävät noin 19 prosenttia koepisteiden vaihtelusta. Selitysaste kertoo tavallaan mallin suorituskyvystä ja se saa arvon nollan ja yhden välillä.
Adjusted R-squared saa aina pienemmän arvon kuin Multiple R-squared sillä siinä huomioidaan mallin monimutkaisuus eli muuttujien määrä. Adjusted R-squared on tarkempi mallin suorituskyvyn mittari kuin Multiple R-squared. Uuden selittävän muuttujan lisääminen malliin kasvattaa todennäköisesti Multiple R-squaredia ja pienentää Adjusted R-squaredia.
par(mfrow = c(2,2))
plot(malli2, which = c(1,2,5))
On piirretty kolme kuvaa, kuva jossa on mallin kaksi sovite ja residuaalit, QQ-plot-kuva jossa on teoreettiset kvantiilit ja standardoidut residuaalit ja kuva jossa on vipuvoima ja standardoidut residuaalit.
Ensimmäisestä kuvasta eli kuvasta jossa on sovite x-akselilla ja residuaalit y-akselilla voidaan tutkia sitä onko virheillä vakio varianssi. Koska havainnot ovat jakautuneet suhteellisen tasaisesti nollasuoran ympärille eivätkä muodosta jotakin muotoa tai mallia, voidaan sanoa että virhetermin varianssi on vakio eikä virheiden suuruus riipu selittävistä muuttujista. Virheiden vakio varianssi on yksi oletus joka usein tehdään kun estimoidaan lineaarinen regressiomalli. Se näyttäisi mallin kaksi tapauksessa olevan voimassa.
QQ-plot-kuvasta eli kuvasta jossa on teoreettiset kvantiilit x-akselilla ja standardoidut residuaalit y-akselilla voidaan tehdä päätelmiä virheiden normaalisuudesta. Koska havainnot ovat QQ-plot-kuvassa suoran päällä, virhetermin voidaan sanoa noudattavan suurin piirtein normaalijakaumaa. Niin usein halutaan olettaa ja nyt mallin kaksi tapauksessa voidaan sanoa että oletus että virheet ovat normaalisti jakaantuneet voidaan pitää voimassa.
Kolmannesta kuvasta, jossa on vipuvoima eli Cookin etäisyys x-akselilla ja standardoidut residuaalit y-akselilla voidaan päätellä vipuvoiman suuruus. Vipuvoima mittaa sitä kuinka paljon yhdellä havainnolla on vaikutusta malliin. Hajontakuva jossa on residuaalit ja vipuvoima voi auttaa tunnistamaan ne havainnot joilla on poikkeuksellisen paljon vaikutusta malliin. Kuvasta nähdään että havainto joka saa suurimman vipuvoiman tai Cookin etäisyyden arvon saa suunnilleen arvon 0.05. Luku on sen verran pieni että vipuvoiman mallissa kaksi voidaan sanoa olevan tavallista tai pientä.
Muita lineaarisen regressiomallin oletuksia ovat virheiden normaalisuuden ja vakion varianssin lisäksi oletukset että virheet eivät korreloi eivätkä riipu selittävistä muuttujista.
alc=read.table("data/alc.txt", sep = ",", header = TRUE)
head(alc)## school sex age address famsize Pstatus Medu Fedu Mjob Fjob
## 1 GP F 18 U GT3 A 4 4 at_home teacher
## 2 GP F 17 U GT3 T 1 1 at_home other
## 3 GP F 15 U LE3 T 1 1 at_home other
## 4 GP F 15 U GT3 T 4 2 health services
## 5 GP F 16 U GT3 T 3 3 other other
## 6 GP M 16 U LE3 T 4 3 services other
## reason nursery internet guardian traveltime studytime failures
## 1 course yes no mother 2 2 0
## 2 course no yes father 1 2 0
## 3 other yes yes mother 1 2 2
## 4 home yes yes mother 1 3 0
## 5 home yes no father 1 2 0
## 6 reputation yes yes mother 1 2 0
## schoolsup famsup paid activities higher romantic famrel freetime goout
## 1 yes no no no yes no 4 3 4
## 2 no yes no no yes no 5 3 3
## 3 yes no yes no yes no 4 3 2
## 4 no yes yes yes yes yes 3 2 2
## 5 no yes yes no yes no 4 3 2
## 6 no yes yes yes yes no 5 4 2
## Dalc Walc health absences G1 G2 G3 alc_use high_use
## 1 1 1 3 5 2 8 8 1.0 FALSE
## 2 1 1 3 3 7 8 8 1.0 FALSE
## 3 2 3 3 8 10 10 11 2.5 TRUE
## 4 1 1 5 1 14 14 14 1.0 FALSE
## 5 1 2 5 2 8 12 12 1.5 FALSE
## 6 1 2 5 8 14 14 14 1.5 FALSEattach(alc)Aineiston muuttujia ovat esimerkiksi alkoholin käyttö viikolla, alkoholin käyttö viikonloppuna, terveys, poissaolot, alkoholin käyttö koko viikon aikana, suuri alkoholinkulutus, opiskeluun käytetty aika, huoltaja, perheen koko, osoite, ikä, sukupuoli ja koulu. Siinä oli vain osa muuttujista. Loput näkyvät tulosteesta.
Tutkin seuraavaksi suuren alkoholinkäytön joka on looginen muuttuja yhteyttä terveyteen, poissaoloihin, opiskeluun käytettyyn aikaan ja sukupuoleen. Hypoteesini on että terveys ja suuri alkoholinkäyttö korreloivat negatiivisesti. Arvioin myös että opiskeluun käytetyllä ajalla ja suurella alkoholinkäytöllä on negatiivinen korrelaatio. Arvioin myös että poissaolot ja suuri alkoholinkäyttö korreloivat positiivisesti. Hypoteesini on että sukupuoli vaikuttaa suureen alkoholinkäyttöön niin että miehissä on enemmän heitä jotka käyttävät paljon alkoholia eli heitä joille high_use-muuttuja saa arvon TRUE.
table(health,high_use)## high_use
## health FALSE TRUE
## 1 35 11
## 2 28 16
## 3 62 19
## 4 49 18
## 5 94 50table(absences,high_use)## high_use
## absences FALSE TRUE
## 0 52 13
## 1 38 13
## 2 42 16
## 3 33 8
## 4 24 12
## 5 16 6
## 6 16 5
## 7 9 3
## 8 14 6
## 9 6 6
## 10 5 2
## 11 2 4
## 12 4 4
## 13 1 1
## 14 1 6
## 16 0 1
## 17 0 1
## 18 1 1
## 19 0 1
## 20 2 0
## 21 1 1
## 26 0 1
## 27 0 1
## 29 0 1
## 44 0 1
## 45 1 0table(studytime,high_use)## high_use
## studytime FALSE TRUE
## 1 58 42
## 2 135 60
## 3 52 8
## 4 23 4table(sex,high_use)## high_use
## sex FALSE TRUE
## F 156 42
## M 112 72barplot(health,high_use)
barplot(absences,high_use)
barplot(studytime,high_use)
boxplot(health,high_use)
boxplot(absences,high_use)
boxplot(studytime,high_use)
boxplot(sex,high_use)
Laatikkokuvista nähdään että hypoteesini terveyden ja suuren alkoholinkäytön ja opiskeluun käytetyn ajan ja suuren alkoholinkäytön välisistä suhteista olivat oikein. Niiden joiden alkoholinkäyttö ei ole suurta terveyden mediaani ja ala- ja yläkvantiilit ovat selvästi korkeammalla tasolla kuin heidän joiden alkoholinkäyttö on suurta vastaavat tunnusluvut. Samoin niiden joiden alkoholinkäyttö ei ole suurta opiskeluun käytetyn ajan mediaani ja ala- ja yläkvantiilit ovat korkeammalla tasolla kuin heidän joiden alkoholinkäyttö on suurta vastaavat tunnusluvut. Alkoholinkäytön kasvaminen siis vaikuttaa terveyteen ja opiskeluun käytettyyn aikaan negatiivisesti. Poissaolojen suhteen arvioni meni pieleen. Heidän joiden alkoholinkäyttö ei ole suurta poissaolojen mediaani ja ylä- ja alakvantiilit ovat korkeammalla tasolla kuin heidän joiden alkoholinkäyttö on suurta vastaavat tunnusluvut. Alkoholinkäytön kasvaminen siis vaikuttaa poissaoloihin negatiivisesti mikä oli päinvastoin kuin mitä arvioin. Ristiintaulukosta nähdään että sellaisia miehiä joiden alkoholinkäyttö on suurta on enemmän kuin sellaisia naisia joiden alkoholinkäyttö on suurta. Sellaisia miehiä joiden alkoholinkäyttö ei ole suurta on vähemmän kuin sellaisia naisia joiden alkoholinkäyttö ei ole suurta. Miehissä on siis oltava suhteellisesti suurempi osuus heitä joiden alkoholinkäyttö on suurta kuin naisissa. Hypoteesini oli sukupuolen ja suuren alkoholinkäytön suhteen osalta oikea.
m = glm(high_use ~ health + absences + studytime + sex, data = alc, family = "binomial")
summary(m)##
## Call:
## glm(formula = high_use ~ health + absences + studytime + sex,
## family = "binomial", data = alc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2074 -0.8367 -0.5983 1.0910 2.1407
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.08010 0.52055 -2.075 0.037995 *
## health 0.04493 0.08631 0.521 0.602675
## absences 0.08881 0.02307 3.850 0.000118 ***
## studytime -0.39817 0.15960 -2.495 0.012603 *
## sexM 0.78234 0.25261 3.097 0.001954 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 465.68 on 381 degrees of freedom
## Residual deviance: 423.00 on 377 degrees of freedom
## AIC: 433
##
## Number of Fisher Scoring iterations: 4coef(m)## (Intercept) health absences studytime sexM
## -1.08009806 0.04492699 0.08881324 -0.39817186 0.78234486Logistisen regressiomallin estimoinnin tulostuksesta nähdään että kun terveys kasvaa yhdellä, logaritmi suuren alkoholinkäytön todennäköisyyden ja yksi miinus suuren alkoholinkäytön todennäköisyyden osamäärästä eli log(p/(1-p)):stä kasvaa 0,04 verran. Samoin jos poissaolot kasvavat yhdellä, log(p/(1-p)) kasvaa noin 0,09 verran. Opiskeluun käytetty aika -muuttujan kerroin tulkitaan vastaavalla tavalla. Sukupuoli on luokiteltu muuttuja ja sen kertoimen tulkinta on erilainen kuin muiden muuttujien kertoimien tulkinta. R on valinnut naiset vertailuryhmäksi ja miehet ovat selittävä muuttuja. Vertailuryhmän eli naisten kerroin on mallin vakiotermi. Miesten kerroin kuvaa eroa vakiotermiin eli naisten kertoimeen. Miesten todellinen kerroin on vakiotermi-miesten kerroin eli noin -1.08-0.78=-0.3. Selittäjistä tilastollisesti merkitsevä kertoimen estimaatti on poissaoloilla, opiskeluun käytetyllä ajalla, miehillä ja vakiotermillä. Se että miesten kertoimen estimaatti on tilastollisesti merkitsevä tarkoittaa että hypoteesi että vakiotermin eli vertailuryhmän kertoimen ja miesten kertoimen erotus on nolla voidaan hylätä yhden prosentin merkitsevyystasolla. Terveyden kertoimen estimaatti ei ole tilastollisesti merkitsevä eli hypoteesia siitä että terveyden kerroin on nolla ei voida hylätä yleisillä merkitsevyystasoilla.
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:GGally':
##
## nasa## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionOR =coef(m) %>% exp
CI=exp(confint(m))## Waiting for profiling to be done...cbind(OR, CI)## OR 2.5 % 97.5 %
## (Intercept) 0.3395622 0.1203826 0.9317189
## health 1.0459515 0.8842249 1.2412355
## absences 1.0928765 1.0468090 1.1460148
## studytime 0.6715466 0.4866129 0.9113044
## sexM 2.1865935 1.3377030 3.6085532Terveyden odds ratio kertoo suuren alkoholinkäytön todennäköisyyden henkilöillä joilla terveys-muuttuja muuttuu yhdellä yksiköllä ja suuren alkoholinkäytön todennäköisyyden henkilöillä jolla terveys-muuttuja ei muutu yhdellä yksiköllä osamäärän. Koska se on suurempaa kuin yksi, terveyden muutos yhdellä yksiköllä liittyy positiivisesti suureen alkoholinkäyttöön. Tämä on eri tulos kuin hypoteesissani arvioin. Myös poissaolojen odds ratio on yli yksi eli poissaolojen muutos yhdellä yksiköllä liittyy positiivisesti suureen alkoholinkäyttöön. Se on hypoteesini mukainen tulos. Opiskeluun käytetyn ajan odds ratio on alle yksi mikä tarkoittaa että opiskeluun käytetyn ajan muutos yhdellä yksiköllä liittyy negatiivisesti suureen alkoholinkäyttöön. Se on arvioni mukainen tulos. Se että on mies näyttää liittyvän positiivisesti suureen alkoholinkäyttöön, sillä miesten odds ratio on yli kaksi. Se on arvioni mukaista.
#Lasketaan malli ilman terveysmuuttujaa jonka kertoimen estimaatti ei ollut tilastollisesti merkitsevä.
m = glm(high_use ~ absences + studytime + sex, data = alc, family =
"binomial")
#Ennustetaan suuren alkoholinkäytön todennäköisyys.
probabilities = predict(m, type = "response")
#Lisätään ennustetut todennäköisyydet aineistoon.
alc = mutate(alc, probability = probabilities)
#Käytetään todennäköisyyksiä suuren alkoholinkäytön ennustamiseen,
alc = mutate(alc, prediction = probabilities>0.5)
#Ristiintaulukoidaan korkea alkoholinkäyttö-muuttuja ennusteiden muuttujasta kanssa.
table(high_use = alc$high_use, prediction = probabilities>0.5)## prediction
## high_use FALSE TRUE
## FALSE 258 10
## TRUE 88 26Ristiintaulukosta nähdään että suuri alkoholinkäyttö-muuttujan arvot poikkeavat siitä logit-mallin avulla tehdyistä ennusteista vähän. 10 kertaa ennuste on suuri alkoholinkäyttö silloin kun muuttuja saa ei suurta alkoholinkäyttöä vastaavan arvon. Samoin 88 kertaa ennuste on ei suurta alkoholinkäyttöä kun muuttuja saa suurta alkoholinkäyttöä vastaavan arvon.
# Haetaan paketit dplyr ja ggplot2.
library(dplyr); library(ggplot2)
# Piirretään hajontakuva korkeasta alkoholin käytöstä ja sen ennustetuista todennäköisyyksistä.
g = ggplot(alc, aes(x = probability, y = high_use, col=prediction))
# Määritellään geom pisteinä ja piirretään uusi kuva.
g+geom_point()
#Taulukoidaan suuri alkoholinkäyttö ja ennusteet siitä.
table(high_use = alc$high_use, prediction = alc$prediction)%>%prop.table()%>%addmargins()## prediction
## high_use FALSE TRUE Sum
## FALSE 0.67539267 0.02617801 0.70157068
## TRUE 0.23036649 0.06806283 0.29842932
## Sum 0.90575916 0.09424084 1.00000000# Määritellään loss function eli keskimääräinen ennustevirhe.
loss_func = function(class, prob) {
n_wrong = abs(class - prob) > 0.5
mean(n_wrong)
}
# Kutsutaan loss funktio jotta voidaan laskea väärien ennusteiden keskimääräinen määrä aineistossa.
loss_func(class = alc$high_use, prob = alc$probability)## [1] 0.2565445Väärin luokiteltuja henkilöitä on noin 26 prosenttia kaikista henkilöistä. Ylemmästä taulukosta nähdään että ennusteen mukaan henkilöitä joilla ei ole suurta alkoholinkäyttöä on noin 91 prosenttia kaikista henkilöistä. Muuttujan arvojen mukaan henkilöitä joilla ei ole suurta alkoholinkäyttöä on noin 70 prosenttia kaikista henkilöistä. Huomataan taas että ennuste eroaa muuttujan arvoista.
``` ***
Clustering and classification
# access the MASS package
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select# load the data
data("Boston")
# explore the dataset
str(Boston)## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...dim(Boston)## [1] 506 14We loaded the Boston data from the MASS package. We explored the structure and the dimensions of the data. From the output of the r command dim we see that the data contains 14 variables and 506 observations. From the output of r command str we see that all variables are numerical and some of them get integer values. From the web page concerning the Boston dataset we saw that the variables are the following
crim
per capita crime rate by town.zn
proportion of residential land zoned for lots over 25,000 sq.ft.indus
proportion of non-retail business acres per town.chas
Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).nox
nitrogen oxides concentration (parts per 10 million).rm
average number of rooms per dwelling.age
proportion of owner-occupied units built prior to 1940.dis
weighted mean of distances to five Boston employment centres.rad
index of accessibility to radial highways.tax
full-value property-tax rate per \$10,000.ptratio
pupil-teacher ratio by town.black
1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town.lstat
lower status of the population (percent).medv
median value of owner-occupied homes in \$1000s.summary(Boston)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08204 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00library(ggplot2)
library(GGally)
pairs(Boston[-1])
p=ggpairs(Boston, mapping = aes(alpha = 0.3), lower = list(combo = wrap("facethist", bins = 20)))
p
# calculate the correlation matrix
cor_matrix<-cor(Boston)
# print the correlation matrix
cor_matrix## crim zn indus chas nox
## crim 1.00000000 -0.20046922 0.40658341 -0.055891582 0.42097171
## zn -0.20046922 1.00000000 -0.53382819 -0.042696719 -0.51660371
## indus 0.40658341 -0.53382819 1.00000000 0.062938027 0.76365145
## chas -0.05589158 -0.04269672 0.06293803 1.000000000 0.09120281
## nox 0.42097171 -0.51660371 0.76365145 0.091202807 1.00000000
## rm -0.21924670 0.31199059 -0.39167585 0.091251225 -0.30218819
## age 0.35273425 -0.56953734 0.64477851 0.086517774 0.73147010
## dis -0.37967009 0.66440822 -0.70802699 -0.099175780 -0.76923011
## rad 0.62550515 -0.31194783 0.59512927 -0.007368241 0.61144056
## tax 0.58276431 -0.31456332 0.72076018 -0.035586518 0.66802320
## ptratio 0.28994558 -0.39167855 0.38324756 -0.121515174 0.18893268
## black -0.38506394 0.17552032 -0.35697654 0.048788485 -0.38005064
## lstat 0.45562148 -0.41299457 0.60379972 -0.053929298 0.59087892
## medv -0.38830461 0.36044534 -0.48372516 0.175260177 -0.42732077
## rm age dis rad tax
## crim -0.21924670 0.35273425 -0.37967009 0.625505145 0.58276431
## zn 0.31199059 -0.56953734 0.66440822 -0.311947826 -0.31456332
## indus -0.39167585 0.64477851 -0.70802699 0.595129275 0.72076018
## chas 0.09125123 0.08651777 -0.09917578 -0.007368241 -0.03558652
## nox -0.30218819 0.73147010 -0.76923011 0.611440563 0.66802320
## rm 1.00000000 -0.24026493 0.20524621 -0.209846668 -0.29204783
## age -0.24026493 1.00000000 -0.74788054 0.456022452 0.50645559
## dis 0.20524621 -0.74788054 1.00000000 -0.494587930 -0.53443158
## rad -0.20984667 0.45602245 -0.49458793 1.000000000 0.91022819
## tax -0.29204783 0.50645559 -0.53443158 0.910228189 1.00000000
## ptratio -0.35550149 0.26151501 -0.23247054 0.464741179 0.46085304
## black 0.12806864 -0.27353398 0.29151167 -0.444412816 -0.44180801
## lstat -0.61380827 0.60233853 -0.49699583 0.488676335 0.54399341
## medv 0.69535995 -0.37695457 0.24992873 -0.381626231 -0.46853593
## ptratio black lstat medv
## crim 0.2899456 -0.38506394 0.4556215 -0.3883046
## zn -0.3916785 0.17552032 -0.4129946 0.3604453
## indus 0.3832476 -0.35697654 0.6037997 -0.4837252
## chas -0.1215152 0.04878848 -0.0539293 0.1752602
## nox 0.1889327 -0.38005064 0.5908789 -0.4273208
## rm -0.3555015 0.12806864 -0.6138083 0.6953599
## age 0.2615150 -0.27353398 0.6023385 -0.3769546
## dis -0.2324705 0.29151167 -0.4969958 0.2499287
## rad 0.4647412 -0.44441282 0.4886763 -0.3816262
## tax 0.4608530 -0.44180801 0.5439934 -0.4685359
## ptratio 1.0000000 -0.17738330 0.3740443 -0.5077867
## black -0.1773833 1.00000000 -0.3660869 0.3334608
## lstat 0.3740443 -0.36608690 1.0000000 -0.7376627
## medv -0.5077867 0.33346082 -0.7376627 1.0000000We showed a graphical overview of the data and summaries of the variables in the data. From the summaries we can see that for example minimum of crime rate in one town is about 0.006, median is about 0.257, mean is about 3.614 and maximum is about 88.976. From the graph we can see that the distributions of one variable against another are rarely normal. Only variable average number of rooms per dwelling´s distribution when the same variable, average number of rooms per dwelling, is on the y-axis looks normal. From correlation matrix we see that for example crime rate in one town correlates a lot with index of accessibility to radial highways and with full-value property-tax rate per 10,000 dollars. Correlation of crime rate with index of accessibility to radial highways is 0.63. Correlation of crime rate and full-value property-tax rate per 10,000 dollars is 0.58.
# center and standardize variables
boston_scaled <- scale(Boston)
# summaries of the scaled variables
summary(boston_scaled)## crim zn indus
## Min. :-0.419367 Min. :-0.48724 Min. :-1.5563
## 1st Qu.:-0.410563 1st Qu.:-0.48724 1st Qu.:-0.8668
## Median :-0.390280 Median :-0.48724 Median :-0.2109
## Mean : 0.000000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.007389 3rd Qu.: 0.04872 3rd Qu.: 1.0150
## Max. : 9.924110 Max. : 3.80047 Max. : 2.4202
## chas nox rm age
## Min. :-0.2723 Min. :-1.4644 Min. :-3.8764 Min. :-2.3331
## 1st Qu.:-0.2723 1st Qu.:-0.9121 1st Qu.:-0.5681 1st Qu.:-0.8366
## Median :-0.2723 Median :-0.1441 Median :-0.1084 Median : 0.3171
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.:-0.2723 3rd Qu.: 0.5981 3rd Qu.: 0.4823 3rd Qu.: 0.9059
## Max. : 3.6648 Max. : 2.7296 Max. : 3.5515 Max. : 1.1164
## dis rad tax ptratio
## Min. :-1.2658 Min. :-0.9819 Min. :-1.3127 Min. :-2.7047
## 1st Qu.:-0.8049 1st Qu.:-0.6373 1st Qu.:-0.7668 1st Qu.:-0.4876
## Median :-0.2790 Median :-0.5225 Median :-0.4642 Median : 0.2746
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.6617 3rd Qu.: 1.6596 3rd Qu.: 1.5294 3rd Qu.: 0.8058
## Max. : 3.9566 Max. : 1.6596 Max. : 1.7964 Max. : 1.6372
## black lstat medv
## Min. :-3.9033 Min. :-1.5296 Min. :-1.9063
## 1st Qu.: 0.2049 1st Qu.:-0.7986 1st Qu.:-0.5989
## Median : 0.3808 Median :-0.1811 Median :-0.1449
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.4332 3rd Qu.: 0.6024 3rd Qu.: 0.2683
## Max. : 0.4406 Max. : 3.5453 Max. : 2.9865# change the object to data frame
boston_scaled=as.data.frame(boston_scaled)We scaled the variables in the Boston dataset. From the summary of the scaled variables we see that for example statistics of crime rate have changed radically. Minimum changed from about 0.006 to about -0.419. Median changed from about 0.257 to about -0.39. Mean changed from about 3.614 to 0 and maximum changed from about 88.976 to about 9.924110.
# save the scaled crim as scaled_crim
scaled_crim <- boston_scaled$crim
# create a quantile vector of crim and print it
bins <- quantile(scaled_crim)
# create a categorical variable 'crime', label-käsky on varmaan väärin
crime <- cut(scaled_crim, breaks = bins, include.lowest = TRUE, label=c("low","med_low","med_high","high"))
# remove original crim from the dataset
boston_scaled <- dplyr::select(boston_scaled, -crim)
# add the new categorical value to scaled data
boston_scaled <- data.frame(boston_scaled, crime)
# number of rows in the Boston dataset
n <- nrow(boston_scaled)
# choose randomly 80% of the rows
ind <- sample(n, size = n * 0.8)
# create train set
train <- boston_scaled[ind,]
# create test set
test <- boston_scaled[-ind,]We have created a categorical variable of the crime rate in the Boston dataset from the scaled crime rate. We used the quantiles as the break points in the categorical variable and dropped the old crime rate variable from the dataset. We also divided the dataset to train and test sets, so that 80% of the data belongs to the train set.
# linear discriminant analysis
lda.fit <- lda(crime ~ ., data = train)
# print the lda.fit object
lda.fit## Call:
## lda(crime ~ ., data = train)
##
## Prior probabilities of groups:
## low med_low med_high high
## 0.2549505 0.2549505 0.2500000 0.2400990
##
## Group means:
## zn indus chas nox rm
## low 0.8673443 -0.8877493 -0.119431971 -0.8473776 0.5187493
## med_low -0.0749120 -0.2837441 -0.004759149 -0.5400482 -0.1770405
## med_high -0.3845049 0.2218593 0.234426408 0.4686402 0.1846551
## high -0.4872402 1.0172187 -0.028797094 1.0657861 -0.4299321
## age dis rad tax ptratio
## low -0.8371977 0.8144318 -0.7020020 -0.7537982 -0.42253114
## med_low -0.3144568 0.3673516 -0.5414397 -0.4614481 -0.06556348
## med_high 0.4552301 -0.4117607 -0.3792103 -0.2635934 -0.35675968
## high 0.8143138 -0.8556025 1.6371072 1.5133254 0.77958792
## black lstat medv
## low 0.37706741 -0.79011874 0.57850729
## med_low 0.34509377 -0.08961643 -0.04389209
## med_high 0.08774016 -0.07077932 0.20915598
## high -0.56847645 0.83822310 -0.65841233
##
## Coefficients of linear discriminants:
## LD1 LD2 LD3
## zn 0.10329131 0.57103194 -0.89545577
## indus -0.04196000 -0.27772378 0.32064043
## chas -0.06970918 -0.03297359 0.12967373
## nox 0.40505960 -0.85246448 -1.11311673
## rm -0.09698941 -0.10402259 -0.21767468
## age 0.33196083 -0.38919579 -0.11651836
## dis -0.08119916 -0.25368979 0.21090130
## rad 3.04446171 0.79463649 -0.11192349
## tax -0.01960874 0.11877117 0.32968885
## ptratio 0.11272780 0.05983056 -0.26340589
## black -0.12325795 0.04555778 0.20014028
## lstat 0.13621171 0.02011664 0.60165781
## medv 0.12561622 -0.30712947 -0.09766526
##
## Proportion of trace:
## LD1 LD2 LD3
## 0.9447 0.0395 0.0158# the function for lda biplot arrows
lda.arrows <- function(x, myscale = 1, arrow_heads = 0.1, color = "orange", tex = 0.75, choices = c(1,2)){
heads <- coef(x)
arrows(x0 = 0, y0 = 0,
x1 = myscale * heads[,choices[1]],
y1 = myscale * heads[,choices[2]], col=color, length = arrow_heads)
text(myscale * heads[,choices], labels = row.names(heads),
cex = tex, col=color, pos=3)
}
# target classes as numeric
classes <- as.numeric(train$crime)
# plot the lda results
plot(lda.fit, dimen = 2, col = classes, pch = classes)
lda.arrows(lda.fit, myscale = 1)
We fitted the linear discriminant analysis that is LDA on the train set. We used the categorical crime rate as the target variable and all the other variables in the dataset as predictor variables. We drew the LDA (bi)plot.
# save the correct classes from test data
correct_classes <- test$crime
# remove the crime variable from test data
test <- dplyr::select(test, -crime)
# predict classes with test data
lda.pred <- predict(lda.fit, newdata = test)
# cross tabulate the results
table(correct = correct_classes, predicted = lda.pred$class)## predicted
## correct low med_low med_high high
## low 17 7 0 0
## med_low 4 17 2 0
## med_high 1 15 9 0
## high 0 0 0 30We saved the crime categories from the test set and then removed the categorical crime variable from the test dataset. Then we predicted the classes with the LDA model on the test data. We also cross tabulated the results from prediction with the crime categories from the test set. We notice that the predicted values of the class variable that are predicted on the test data using the results from the LDA on the train data differ partly from the observed values of the class variable in the test data. Anyhow the predictions are same as the observed values more frequently that the predictions differ from observed values. Prediction is hence quite good.
# load the data
data("Boston")
# center and standardize variables
boston_scaled <- scale(Boston)
# change the object to data frame
boston_scaled=as.data.frame(boston_scaled)
# euclidean distance matrix
dist_eu <- dist(boston_scaled)
# look at the summary of the distances
summary(dist_eu)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1343 3.4620 4.8240 4.9110 6.1860 14.4000# manhattan distance matrix
dist_man <- dist(boston_scaled, method="manhattan")
# look at the summary of the distances
summary(dist_man)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2662 8.4830 12.6100 13.5500 17.7600 48.8600# k-means clustering
km <-kmeans(dist_eu, centers = 15)
# plot the Boston dataset with clusters
pairs(boston_scaled, col = km$cluster)
set.seed(123)
# determine the number of clusters
k_max <- 15
# calculate the total within sum of squares
twcss <- sapply(1:k_max, function(k){kmeans(dist_eu, k)$tot.withinss})
# visualize the results
plot(1:k_max, twcss, type='b')
# k-means clustering
km <-kmeans(dist_eu, centers = 2)
# plot the Boston dataset with clusters
pairs(boston_scaled, col = km$cluster)
We reloaded the Boston dataset and standardized the dataset. We calculated the distances between the observations using the euclidian and manhattan distance matrixes. We ran k-means algorithm on the dataset and visualized the clusters with pairs-function where clusters are separated with colors. Then we investigated what is the optimal number of clusters is with total within sum of squares. We drew a picture where the number of clusters in on the x-axis and quantity of total within sum of squares on the y-axis. From the picture we saw that total sum of squares gets its highest value at one. We thought one is too small number of clusters and chose two for the number of clusters. We ran the k-means algorithm again and chose two as the number of centers. We visualized the clusters again with pairs-function.
Super-bonus-exercise
model_predictors <- dplyr::select(train, -crime)
# check the dimensions
dim(model_predictors)## [1] 404 13dim(lda.fit$scaling)## [1] 13 3# matrix multiplication
matrix_product <- as.matrix(model_predictors) %*% lda.fit$scaling
matrix_product <- as.data.frame(matrix_product)
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutplot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers')plot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers', color=train$crime)